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In the context of the exact factorization of the electron-nuclear wave function, the coupling be¬ 
tween electrons and nuclei beyond the adiabatic regime is encoded (i) in the time-dependent vector 
and scalar potentials and (ii) in the electron-nuclear coupling operator. The former appear in the 
Schrodinger-like equation that drives the evolution of fhe nuclear degrees of freedom, whereas the 
latter is responsible for inducing non-adiabatic effects in the electronic evolution equation. As we 
have devoted previous studies to the analysis of the vector and scalar potentials, in this paper we fo¬ 
cus on the properties of the electron-nuclear coupling operator, with the aim of describing a numerical 
procedure fo approximate it within a semiclassical treatment of the nuclear dynamics. 


I. INTRODUCTION 

Modelling the d 5 rnamical coupling of electrons and 
nuclei beyond the Bom-Oppenheimer (BO), or adia¬ 
batic, regime is currently among the most challenging 
problems in the fields of Theoretical Chemistry and 
Condensed Matter Physics. Within the BO framework, 
molecular systems are visualized as a set of nuclei mov¬ 
ing on a single potential energy surface that represents 
the effect of the electrons in a given eigenstate. Many 
interesting phenomena, however, such as vision fl] |2l, 
charge separation in organic photovoltaic materials ||3l 
m or Joule heating in molecular junctions |]5l |6l, occur 
in non-adiabatic conditions. In these situations, solv¬ 
ing exactly the time-dependent Schrodinger equation 
(TDSE) for the coupled system of electrons and nuclear 
is not feasible, as the computational cost scales expo¬ 
nentially with the number of degrees of freedom. How¬ 
ever, since the full quantum treatment requires to repre¬ 
sent the problem in terms of adiabatic states and tran¬ 
sitions among them in regions of strong non-adiabatic 
coupling, wave packet propagation techniques have 
been developed, retaining the quantum character of the 
nuclear dynamics 1171 120] . And these techniques are 
presently the state of the art in quantum d 5 mamics com¬ 
putational methods, proving the benchmark for approx¬ 
imate methods. In fact, for large systems, the dimen¬ 
sionality of the problem does not allow to employ a 
quantum mechanical description, thus the only feasible 
approach is to combine a classical description for the nu¬ 
clei (or ions) with a quantum treatment of a few other es¬ 
sential degrees of freedom, e.g. electrons or protons . In 
this context, the question of how to model the coupling 
between the quantum and classical subsystems still re¬ 
mains open, despite the fact that several schemes EOl - 
|3^ have been proposed in the literature trying to settle 
this issue. 

In recent work, we have addressed this problem in 
the context of the exact factorization of the electron- 
nuclear wave function Il40l 111] . In such a treatment of 
quantum d 5 mamics, the solution of the TDSE is written 


as a single product of a nuclear wave function and an 
electronic factor, that parametrically depends on the nu¬ 
clear configuration. Several advantages of this reformu¬ 
lation have been pointed out. Eirst of all, it has been 
shown Il42l |43] that the nuclear wave function evolves 
according to a modified TDSE where a time-dependent 
vector potential and a time-dependent scalar potential 
represent the effect of the electrons on the nuclei, be¬ 
yond the adiabatic regime. When a classical treatment 
of the nuclear degrees of freedom is introduced, the cou¬ 
pling to the electrons is exactly represented by the force 
determined from the gradient of the time-dependent po¬ 
tentials Il43ll44] . On the other hand, the electronic factor 
evolves according to a (less standard) evolution equa¬ 
tion, coupled to the nuclear TDSE, where the effect of 
the nuclei is represented by an electron-nuclear coupling 
operator explicitly depending on the nuclear wave func¬ 
tion. 

The analysis of the time-dependent potentials and of 
the classical nuclear force has been the subject of pre¬ 
vious work. We have been able to analyze a simple 
model system to pinpoint some relevant features Il42l - 
|44] of the potentials that should be accounted for when 
developing approximations. Based on these observa¬ 
tions, starting from the exact formulation, we have pro¬ 
posed Il4^ |46] a novel mixed quantum-classical algo¬ 
rithm to solve the coupled electronic and nuclear evo¬ 
lution equations in a fully approximated way. The clas¬ 
sical limit is considered as the lowest-order, in a h- 
expansion, of the nuclear wave function in the complex- 
phase representation 147] . 

In the present paper, the focus is directed towards 
the analysis of the electron-nuclear coupling operator 
that, in the electronic equation, mediates the coupling 
to the nuclei and that depends explicitly on (the gra¬ 
dient of) the nuclear wave function. It is fundamen¬ 
tal to be able to correctly approximate such term, as it 
is responsible to induce electronic non-adiabatic transi¬ 
tions Il4^l46l and decoherence 148] . When a classical de¬ 
scription of the nuclei is adopted, the concept of wave 
function is somehow lost and problems arise when ap- 


proximating this operator. If a distribution of frajecfo- 
ries 144] is used fo mimic fhe evolution of fhe nuclear 
wave function, ifs modulus and phase carmof be smoofh 
functions of space. The numerical error fhus infroduced 
affecfs fhe calculations, buf if can be cured if refined ap¬ 
proximations are considered. The goal of fhis paper is 
to describe a procedure to avoid the above issue and to 
test its efficiency against exact calculations for a simple 
model system. Therefore, we propose here (i) to employ 
a representation of the nuclear density in terms of evolv¬ 
ing frozen gaussians (FGs) l49l , rather than trajectories, 
following the scheme presented in Ref. l44] and (ii) to 
estimate the phase of the nuclear wave function adopt¬ 
ing such a FGs picture, based on a simplified form of the 
semiclassical Herman-Kluk l50][53] propagation scheme 
within the initial value representation (IVR) theory l54l - 
[57l . The model system for non-adiabatic charge trans¬ 
fer of Ref. f58l allows for an exact numerical solution of 
the full quantum mechanical problem, thus providing 
a benchmark to any approximation that will be consid¬ 
ered. Starting from these results, we will compute the 
exact time-dependent scalar potential, also referred to 
as time-dependent potential energy surface (TDPES), in 
a gauge where the vector potential can be set to zero. 
The effect of the electrons is then fully accounted for by 
this TDPES, that is adopted to evolve the EGs. Since the 
electronic part of the problem is solved exactly, the only 
source of error will be in this semiclassical approxima¬ 
tion. It is worth stressing that this procedure does not 
result in the development of a new algorithm, but is a 
test of the performance of the EGs approximation of the 
nuclear motion. 


we 


The paper is organized as follows. In Section Iffl 
briefly recall the factorization formalism and we focus 
on the analysis of the electron-nuclear coupling (ENG) 


operator in the electronic evolution equation. Section III 


is devoted to a discussion on the apronximations em¬ 
ployed in the calculations. Numerical results are pre¬ 
sented in Section |IV1 by comparing situations with dif¬ 


ferent non-adiabatic coupling strengths and testing the 
approximations employed to evaluate the ENG term. 
Gonclusions are stated in Section Ivl 


II. THE EXACT EACTORIZATION ERAMEWORK 


can be written as the product 




( 2 ) 


of the nuclear wave function, x(Ej the electronic 

wave function, 4’R(r, t), which parametrically depends 
of the nuclear configuration Il59ll60] . Throughout the pa¬ 
per the symbols r, R indicate the coordinates of the Ng 
electrons and 7V„ nuclei, respectively. Eq. 0 is unique 
under the partial normalization condition (PNC) 
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up to within a gauge-like phase transformation. The 
evolution equations for <l>R(r, t) and x(Ej t)/ 

(-ffeZ - e(E, t)) t) = t) (4) 

HnXiR, t) = t), (5) 

are derived by applying Erenkel's action principle IIM] 
l62! with respect to the two wave functions and are ex¬ 
actly equivalent lEOlIlT] to the TDSEffl. Eqs. Q and ||^ 
are obtained by imposing the PNC 1631164] by means of 
Lagrange multipliers. 

The electronic equation Q contains the electronic 
Hamiltonian 


Hgi =Hbo (t, R) + [‘Fe, x] , (6) 

which is the sum of the BO Hamiltonian and the ENG 
operator {7 ““p[$r, x]. 
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-|-Ai,(R,f) j - A„(R, t)) 


In Eq. (|^, e(R, t) is the TDPES, defined as 


e(R,t) = (4-R(t) 


Hg! ihdf. 




( 8 ) 


and e(R, t), along with the vector potential 
A(R,t), 


In the absence of an external field, the non-relativistic 
Hamiltonian H = + Hbo describes a system of in¬ 

teracting nuclei and electrons. Here, denotes the nu¬ 
clear kinetic energy and Hbo{^-, E) = Tg{r) -I- 14.n(r, E) 
is the BO Hamiltonian, containing the electronic kinetic 
energy Tg(r) and all interactions l7e_„(r,R). As recently 
proven 14011411 , the full wave function, T'(r, R, t), solu¬ 
tion of the TDSE 


tT4'(r, R, t) = t/idt'I'(r, R, t), (1) 


A(R,t) = ($R(t) -iW^$R(t) 


(9) 


mediate the coupling between electrons and nuclei in a 
formally exact way. Here, the symbol (■ | • )r stands for 
an integration over electronic coordinates. 

The nuclear evolution is generated by the Hamilto¬ 
nian 


Hn{R,t)=Y. 


[-zSV, + A,(R,t)]" 
2A4 


e(R,t), (10) 
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according to the TDSE ||^. 

The TDPES and the vector potential are uniquely de¬ 
termined up to within gauge-like transformations Il40l 
|4TI . The uniqueness can be straightforwardly proven by 
following fhe sfeps of fhe currenf densify version 1651 
of fhe Runge-Gross fheorem USSli . In fhis paper, as a 
choice of gauge, we infroduce fhe addifional consfrainf 
Ai, (R, t) = 0 (see Ref. lilll for a defailed discussion on 
how this condition can be imposed) l67l . 

As discussed in the introduction, we study the prop¬ 
erties of fhe ENC ferm —i^vxlx that in the expres¬ 
sion of fhe operafor 17g°“P[<l>R, y] explicifly depends on 
fhe nuclear wave funcfion. This analysis is based on 
fhe inferesf in developing a procedure fo approximate 
it when a classical or semiclassical treatment of fhe nu¬ 
clear mofion is adopfed. Eor insfance, in fhe classical 
limif, we have derived Il45ll46l ifs expression in ferms of 
fhe nuclear momenfum. This has been done by wrifing 
fhe nuclear wave funcfion as x(R, f) = exp [|5(R, f)], 
wifh 5(R, f) a complex funcfion 1471 . If we now sup¬ 
pose 14511461 fhat fhis funcfion can be expanded as an 
asympfofic series in powers of h, namely S (R, f) = 
Sa (R, t), fhe ENC ferm becomes 


-ffiV„x(R,f) 

x(E,f) 


= V„S'o(R,f), 


( 11 ) 


af fhe lowesf-order in H. On fhe righf-hand-side (RHS), 
fhe funcfion VyS'o(R, t) is fhe classical nuclear momen¬ 
fum evaluafed along fhe frajecfory, since l45ll46l Sq saf- 
isfies a Hamilfon-Jacobi equafion wifh Hamilfonian 


Hr 


^ [y.So {R,t) + A,{R,t)y 


2 M„ 


+ e(R,f). (12) 


The vecfor pofenfial appears in fhe above expression of 
fhe classical Hamilfonian because fhis resulf has general 
validity, not only in the gauge adopted in the following 
calculations. 

Alternatively, if fhe nuclear wave funcfion is written 
in ferms of ifs modulus and phase, x = the (ex- 

acf) expression of fhe ENC ferm becomes 


-fhV„x(R,f) 


= V^5'(R,f)-hf 


|x(R,f) 

|x(E, t)| 


(13) 


If is clear af fhis poinf fhaf a good esfimafe of fhe ENC 
ferm is only possible when bofh fhe modulus and fhe 
phase of fhe nuclear wave funcfion are correcfly de¬ 
scribed. A classical freafmenf, as in Eq. | [TT) , only pro¬ 
vides an approximation fo fhe real parf of fhe ENC 
ferm, while fhe informafion abouf fhe imaginary parf is 
losf IMl . 

In fhe following, we will infroduce a EC-based ap¬ 
proach fo defermine an approximation to Eq. (T^ . 
EGs ||49^ evolving on the exact TDPES are used to re¬ 
construct the nuclear density, thus allowing to calculate 


the second term on the RHS, as |x| is a smooth function 
of fhe nuclear coordinafes. The phase informafion is in- 
sfead encoded in fhe classical action accumulafed over 
time and associated to each EG, as we will show in Sec¬ 
tion |nl] 

Henceforth, we will drop the bold-double underlined 
notation for elecfronic and nuclear positions as we will 
deal wifh one-dimensional (ID) quantifies. 


III. SEMICLASSICAL APPROXIMATION 


A. Nuclear density 


According to the procedure presented in Ref. 11441 , a 
sef of independenf classical frajecfories evolving on the 
exact TDPES are able to reproduce the nuclear density 
in almost perfect agreement with quantum results. As 
pointed out in the Introduction, however, constructing a 
histogram from fhe disfribufion of fhe frajecfories does 
nof allow fo compufe fhe second ferm on fhe RHS of 
Eq. 1131, fhaf involves fhe gradient of |x|, wifhouf a large 
numerical error. The reason is fhaf fhe "classical" den¬ 
sify is nof a smoofh funcfion of space. The solution pro¬ 
posed here is fo improve fhe previous approximafion of 
nuclear d 3 mamics, by propagating fhe mean positions 
and momenfa of a sef of EGs on fhe exacf TDPES, rafher 
fhan classical frajecfories. 

Given a sef of Nt^aj initial positions and momenfa, 
Ro, Pq, sampled as described in Secfion complex 
gaussians, also referred fo as coherent states, are con- 
sfrucfed as 


giR;Ri,oPLon) = (I) 

(14) 

wifh widfh 7 fo be defermined below. Each EG is also 
associafed a "weigh!". 


= y dRg* {R; i?z,oR,o, 7 ) Xo(f?), (15) 

corresponding fo fhe projecfion of fhe initial nuclear 
wave funcfion xo(f?) on fhe EGs. The nuclear densify 
af each time is fhen obfained as 

Ntraj 

\x{R,t)? - XI , (16) 


where Ri{t),Pi{t) are fhe time-evolved positions and 
momenfa of Ri^, Pi^. In comparison fo fhe purely clas¬ 
sical approximafion, Eq. (T^ allows nof only fo repro¬ 
duce fhe nuclear densify in very good agreemenf wifh 
quanfum resulfs, as will be shown in Secfion IV buf 
also fo calculafe (analyfically) fhe gradient of fhe nuclear 
densify (or of fhe modulus, as if appears in Eq. (13)). 

If is imporfanf fo notice fhaf Eq. (T^ is an approxi¬ 
mafion fo fhe nuclear densify when the nuclear wave 
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function is represented as a superposition of coherent 
states. In fact, coherent states form an overcomplete ba¬ 
sis and, in writing Eq. 1 16 1 , we neglect the overlaps of 
coherent states. The reason for this further approxima¬ 
tion is related to the choice of the initial set of positions 
and momenta, the mean positions and mean momenta 
of the FGs. On one hand, we want to maintain the same 
choice of initial conditions done for the classical prop¬ 
agation (see Section |IV[ |, in order to be able to directly 
compare classical arid EG results. On the other hand, 
we want to obtain an initial nuclear density as close as 
possible to the exact density. We have thus computed 
the root mean square deviation (RMSD) for the nuclear 
density when either (1) Eq. {1^ is employed or (2) the 
full expression is consideredTi.e. with overlap terms. 
The best agreement achieved in case (1) is for 7 = 7.0 a, 


-2 


(used in Section [TV), with RMSD = 0.020,0.017, 0.006 for 
= 2000, 5000, 10000, respectively, whereas in case 


Ni 


traj 


(2) is for 7 = 3.0 a^^ with RMSD = 0.037, 0.024 0.0013. 
It is evident that a better agreement at the initial time 
is obtained for case ( 1 ), namely when the approxima¬ 
tion in Eq. | [T^ is used, for all the values of Ntraj- Once 
the parameters defining the coherent states are selected, 
namely Ri, Pi, 7, the weights wi associated to each EG 
are automatically determined by Eq. •ED and kept con¬ 
stant throughout the propagation. 


B. Nuclear phase 


The phase of the nuclear wave function will be deter¬ 
mined according to 


S{R, t) ~ arctan 


'^[xsc{R,t)] 
[xsciR,t)] ’ 


(17) 


where xsc{RR) is a semiclassical approximation to 
the exact x(i?, t). Following the Herman-Kluk proce¬ 
dure ll50H53l to approximate the quantum propagator, 
the expression of the time-evolved wave function at 
time t is 


use a FGs approximation rather than a rigorous semi¬ 
classical approximation. 

In the procedure employed here, the semiclassical nu¬ 
clear wave function is estimated as a sum over trajecto¬ 
ries of time-evolved coherent states, namely 

^traj —Si 

XFG{R,t)= wi^^g{R;Ri{t),Pi{t),x), (19) 

where the /-th classical action is calculated as 

Here, the (exact) TDPES is evaluated, at time r, at the 
classical position Ri (t) and its expression is obtained by 
solving the TDSE for the full wave function T* and by 
directly calculating Eq. ||^ once the factorization ||^ is 
applied. 


IV. NUMERICAL RESULTS 


The expression of the TDPES according to Eq. ||^ is 
determined by calculating the electronic wave function 
$/j(r, t) from the full wave function 'I'(r, R, t), which is 
known at all times by solving the TDSE Q for the model 
Hamiltonian 


Hir,R) 


1 52 1 52 1 

2~^~ ^ \^-R 


1 

U + r\ 

( 21 ) 


e.t(^) erf(h7l) 

|i?-r| k-#! 



This system has been introduce by Shin and Metiu II58II 
as a protot 5 q)e for non-adiabatic charge transfer. The 
system is ID and consists of three ions and a single elec¬ 
tron, as depicted in Fig. Two ions are fixed at a dis- 


xsMR.t) = / 

{R\RtPt,x) {RoPo,l\xo) ■ 

(18) 

Here, (i?oTo,7lxo) denotes the projection of the initial 
nuclear wave function on the coherent states, similarly 
to Eq. | [T^ , while {R\RtPt,x) is an alternative expres¬ 
sion forme coherent states (see Eq. |[T4)). Rt,Pt are 
the (classically) evolved positions and momenta corre¬ 
sponding to the initial conditions Ro,Po- S{Ro, Po',t) 
is the classical action accumulated up to time t along 
the trajectory whose initial conditions are Ro,Po. The 
Herman-Kluk pre-factor is indicated here with the sym¬ 
bol Ct{Ro, Po) II 50 H 53 I , but it will be set equal to unity 
throughout the calculations. Therefore, the symbol 
Xsc{R,t) will be replaced by XFG{R,'t)r since we will 


R 


, . ion 0 electron . 

iixed ion iixea ion 


L 

FIG. 1. Schematic representation of the model system de¬ 
scribed by the Hamiltonian (21) . 

tance of T = 19.0 oq, the third ion and the electron are 
free to move in ID along the line joining the two fixed 
ions. Here, the symbols r and R are the coordinates 
of the electron and the movable ion measured from the 
center of the two fixed ions. The ionic mass is chosen as 
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M = 1836, the proton mass, whereas the other parame¬ 
ters are tuned in order to make the system essentially 
a two-electronic-state model. We present here the re¬ 
sults obtained by choosing two different sets of param- 
efers, producing sfrong and weak non-adiabafic cou¬ 
plings, between fhe firsf, e^Q, and fhe second BOPES, 
( 2 ) 

e^Q, around fhe avoided crossing af Rac = —1-90 ao- 
The values of fhe paramefers in fhe Hamilfonian | [2T) 
are: Rf = 5.0 ao, Ri = 3.1 ao and R^ = 4.0 ao, for 
fhe sfrong coupling case; Rf = 3.8 ao, Ri = 2.9 ag and 
Rr = 5.5 ao, for fhe weak coupling case. The BO surfaces 
are shown in Fig.|^(upper panels). 


strong coupling weak coupling 



R (ao) R (ao) 



0 10 20 30 0 10 20 30 

time (fs) time (fs) 

FIG. 2. Upper panels: Lowest four BO surfaces, as functions of 
the nuclear coordinate, for strong (left) and weak (right) non- 
adiabatic coupling strengths. The first (red line) and second 
(green line) surfaces correspond to the adiabatic states that are 
populated during the dynamics, whereas the third and fourth 
surfaces (dashed black lines) are shown for reference. The 
squared modulus (reduced by ten times and rigidly shifted 
in order to superimpose it on the energy curves) of the ini¬ 
tial nuclear wave packet is also shown (thin black line). Lower 
panels: Populations of the two lowest adiabatic states (pi and 
P 2 ) as functions of time. The arrows represent the time-steps 
shown in the following figures. 


We sfudy fhe time evolution of fhis sysfem by choos¬ 
ing fhe initial wave function as fhe producf of a real¬ 
valued normalized Gaussian wave packef, cenfered af 
Rc — —4.0 ao wifh variance tr = l/-\/2.85 ao (fhin 
black line in Fig. 0 upper panels), and fhe second BO 
elecfronic sfafe. To calculafe fhe TDPES, we firsf solve 
fhe TDSE iQ for fhe complefe sysfem, wifh Hamilfo¬ 
nian | [2T) , and obfain fhe full wave function, T'(r, R, t). 
This isdone by numerical infegrafion of fhe TDSE us¬ 
ing fhe splif-operafor-fechnique 16911 , wifh fime-sfep of 
2.4 X 10“^ fs (or 0.1 a.u.). 

The mean positions and momenfa of fhe FGs evolve 
along classical frajecfories, generafed according fo 


Hamilfon's equafions 


Ri{t) = Pi{t)/M- Pi(t) = -dneiR,t) (22) 

infegrafed by using fhe velocify-Verlef algorifhm wifh 
a fime-sfep of 12.0 x 10“^ fs (or 0.5 a.u.). The initial 
conditions are sampled from fhe Wigner phase-space 
disfribufion corresponding fo fhe inifial nuclear wave 
function. The inifial coherenf sfafes used in Eq. | [T5| | 
are fhus consfrucfed fo defermine fhe (complex-valueH) 
weighfs wi. Numerical resulfs have been obfained for 
sefs of Ntraj = 2000,5000,10000 frajecfories, wifh fhe 
value iTOl 7 = 7.0 ajj"^ for fhe widfh of fhe FGs. The 
coherenf sfafes used fo represenf fhe nuclear wave func¬ 
tion form an overcomplefe basis, fherefore fhe sum in 
Eq. ( [Th) has fo be fruncafed. In order fo choose fhe ad- 
equafe number of basis funcfions fo include in fhe sum, 
energy conservafion has been fesfed and confirmed for 
all values of Ntraj- The resulfs will be presenfed only for 
fhe case Ntraj = 5000. We have compufed fhe RMSD be- 
fween fhe exacf nuclear densify and ifs approximation 
in Eq. ( [Th) for a sef of values of 7 and we have chosen 
fhe value of fhis paramefer for which fhe RMSD is mini¬ 
mum (see also fhe discussion af fhe end of Secfion III A} , 
buf clearly ofher values can be selecfed. 

We will confirm below fhaf fhe semiclassical EG 
scheme, adapfed fo fhe facforizafion approach, is an 
accurafe and efficienf way fo approximafe fhe nuclear 
wave function. Before presenfing fhe EG resulfs, lef us 
firsf show differenf snapshofs faken along fhe d 5 mam- 
ics, showing fhe TDPES. Ifs feafure have been exfen- 
sively discussed in previous work 14211441 . buf for fhe 
sake of complefeness, we reporf here a few configura- 
fions. Moreover, we would like fo underline fhaf here re¬ 
sulfs for differenf non-adiabafic coupling sfrengfhs will 
be presenfed, in order fo fesf fhe efficiency of fhe mefhod 
also in fhe weak coupling regime. 

Fig. shows fhe gauge-invarianf (GI) and gauge- 
dependenf (GD) componenfs of fhe TDPES, namely 
fhe fwo ferms fhaf can be idenfified in Eq. as 
eG/(S,f) = (4’R(f)|fTez|4>R(t))r and eon {^,t) = 

(4>R(f)| — ifidt\^u{t))r- The snapshofs shown in Fig.j^ 
are faken along fhe evolufion af fhe times indicafed by 
fhe arrows in Fig. (lower panels) and fhe fwo lowesf 
adiabafic surfaces are shown for reference. As previ¬ 
ously discussed l42ll44l , before fhe spliffing of fhe nu¬ 
clear wave packef, fhe GI parf of fhe TDPES is diabafic, 
whereas fhe GD parf is consfanf, and affer fhe splif¬ 
fing fhe GI componenf develops sfeps fhaf bridge be- 
fween differenf adiabafic surfaces, whereas fhe GD parf 
is piecewise consfanf. The TDPES is only calculafed in 
fhe regions where fhe nuclear densify is (numerically) 
nof zero. The lack of reliable information beyond fhe re¬ 
gions shown in fhe figures is nof an issue when classical 
frajecfories or FGs are employed fo mimic fhe nuclear 
densify, as fhe regions where fhe exacf densify is expo¬ 
nentially small are nof, or are poorly, sampled. 

If is evidenf from Fig.|^fhaf we will discuss resulfs for 
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FIG. 3. TDPES at different time-steps, as indicated by the ar¬ 
rows in Fig.|^ for sfrong (left) and weak (right) non-adiabatic 
coupling strengths. The two lowest BO surfaces (black) are 
plotfed for reference. The GI part of fhe TDPES (red) presents 
the well-studied |42] IHl dynamical steps that bridge piece- 
wise adiabatic shapes potential energy surface, whereas the 
GI part l44l (cyan) is piecewise constant. The nuclear density 
(blue) is shown for reference. 


short d 5 mamics, limited to the first half of the oscillation 
period of the nuclear wave packet in the potential well. 
Interesting d 5 mamics may arise at later times, when for 
instance the nuclear wave packet crosses a second time 
the non-adiabatic coupling region. However, here we 
focus on the initial non-adiabatic event and test how the 
FG approximation capture this process. Due to the fact 
that the analysis reported below is the first attempt to 
incorporate semiclassically nuclear quantum effects in 
the exact factorization formalism, we study a simple sit¬ 
uation that nonetheless captures the main features of 
a non-adiabatic event. Also, situations where, for in¬ 
stance, reproducing tunneling d 5 mamics might repre¬ 
sent a problem for classically evolving FGs ll7TI - (74l will 
be the subject of further study and not addressed here. 

A first set of results is shown in Fig. where we 
compare the nuclear density, |x(i?, t)p, from three dif¬ 
ferent calculations: quantum (cyan), employing the full 
electron-nuclear wave function; classical (green), where 
the histogram is constructed from the distribution of 
classical trajectories evolving on the TDPES according 
to Eqs. I p^ (as in Ref. l44l ): EG (red), with the nu¬ 
clear density given in Eq. ( [T^ . As expected from pre¬ 
vious calculations 1441 , the use of classical trajectories 
seems to be enough accurate to reproduce the nuclear 
density. It is however important to stress again, that 
these results are not obtained by solving a fully approx¬ 
imate form of the coupled electronic and nuclear equa- 
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FIG. 4. Nuclear density for strong (left) and weak (right) non- 
adiabatic coupling strengths. Exact results (cyan) are com¬ 
pared with the semiclassical density (red), expressed as sum 
of EGs, and with the histogram (green) constructed from the 
distribution of classical positions along the trajectories. 


tions Q and ||5l. They only represent a benchmark for 
any quantum-aassical algorithm, since the effect of the 
electrons, via the TDPES, is treated exactly. The semi¬ 
classical density, constructed as the weighted sum of 
EGs given in Eq. (T^ , is also accurate. In comparison to 
the classical histogram, the gain here is the smoothness 
of the density, not achievable with purely classical tra¬ 
jectories. This feature is extremely important for the cal¬ 
culation of the ENG term containing the gradient of the 
nuclear density, via the term —iKS/u\x\/\x\ Eq. | |T^ . 

The second important characteristic of the semiclassi¬ 
cal approach is that each EG contributes a phase factor 
to the full nuclear wave function, thus allowing to deter¬ 
mine also the first term on the RHS of Eq. < [T^ . We show 
this term, i.e. V^S{R, t), in Eig.[^ comparing once again 
exact results with the corresponding EG and classical 
approximations. The semiclassical value of \7^S{R,t) 
is determined as the gradient of the phase in Eq. | [T7| , 
whereas the function Vi,S{R, t) is classically interpreted 
as the nuclear momentum evaluated along each trajec¬ 
tory. While the agreement between quantum and EG re¬ 
sults is remarkable, the phase-space points correspond¬ 
ing to the classical trajectories do not allow to recon¬ 
struct a smooth function of R. Even if the number of 
classical trajectories is increased, the phase-space points 
are too "noisy" to allow for reconstructing a smooth 
function. A smoothing algorithm should then be em¬ 
ployed, but the numerical efficiency of the whole proce¬ 
dure might become questionable. 

EigJ^shows the imaginary part of the ENG term from 
Eq. | jl3^ at different time-steps during the d 3 mamics, as 
in previous figures. The exact results (cyan) are com- 
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FIG. 5. Real part of the ENC term for strong (left) and weak 
(right) non-adiabatic coupling strengths. The nuclear density 
(black) is shown for reference. Exact results (cyan) are com¬ 
pared with semiclassical calculations (red) and with the classi¬ 
cal phase-space points (green). 


strong coupling weak coupling 


10 n-1-1- 

t = 4.84 fs FG 
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FIG. 6. Imaginary part of the ENC term for strong (left) and 
weak (right) non-adiabatic coupling strengths. The color code 
is the same as in Fig.|^ 


pared only with the approximation (red) based on the 
semiclassical propagation of FGs. Even if we employ a 


simplified form of the Herman-Kluk propagator, where 
the pre-factor is set to 1, semiclassical results are in sat¬ 
isfactory good agreement with exact results. 

One could consider improving the approach devel¬ 
oped here by, for instance, explicitly computing the 
Herman-Kluk pre-factor, at the expenses of increasing 
the computational cost. 


V. CONCLUSION 

We have reported our first semiclassical procedure 
adapted to the formalism of the exact factorization of the 
electron-nuclear wave function lllOl IITI . The approach 
has been used to estimate the ENC term that explic¬ 
itly depends on the nuclear wave function (gradient of 
its modulus and phase) in the electronic equation. In 
previous work 14^ l46ll on the development of a mixed 
quantum-classical algorithm in the context of the exact 
factorization, such term has been treated fully classically 
and identified as the nuclear momentum. However, we 
observed that, despite the fact that this approximation is 
widely used in the literature Il23ll26ll27l , correction terms 
naturally arise in the factorization framework. The ex¬ 
tent and role of the corrections may depend on the pro¬ 
cess that we intend to study, therefore we have shown 
in the present paper how to evaluate such corrections 
based on the semiclassical propagation of EGs, for a sim¬ 
ple model case of non-adiabatic charge transfer. The 
main gain in using EGs, rather than a purely classical ap¬ 
proach, lies in the possibility of obtaining smooth func¬ 
tions whose gradients can be easily determined without 
introducing large numerical errors. 

The semiclassical results shown in the paper have 
been obtained by evolving EGs on the exact TDPES, 
that is known for the simple model system studied here 
as the outcome of exact calculations based on the nu¬ 
merical solution of the full TDSE. This procedure is a 
test, since the only approximation is the semiclassical 
treatment of the nuclear d 3 mamics, whereas the elec¬ 
trons are treated exactly, via the information encoded 
in the TDPES. Moreover, this study provides a bench¬ 
mark for future development, aiming at improving the 
initial, and lowest-order, mixed quantum-classical algo¬ 
rithm derived from the factorization Il4^l46l . The proce¬ 
dure described here can be easily implemented in an al¬ 
gorithm, as we will present elsewhere Il48ll , resulting in a 
novel mixed quantum-semiclassical scheme for solving 
coupled electron-nuclear d 5 mamics. 
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